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ABSTRACT 

Beyond convergence studies and comparison of different codes, there are essentially no 
controls on the accuracy in the non-linear regime of cosmological N body simulations, 
even in the dissipationless limit. We propose and explore here a simple test which has 
not been previously employed: when cosmological codes are used to simulate an isolated 
overdensity, they should reproduce, in physical coordinates, those obtained in open 
boundary conditions without expansion. In particular, the desired coUisionless nature 
of the simulations can be probed by testing for stability in physical coordinates of 
virialized equilibria. We investigate and illustrate the test using a suite of simulations 
in an Einstein de Sitter cosmology from initial conditions which rapidly settle to virial 
equilibrium. We find that the criterion of stable clustering allows one to determine, for 
given particle number N in the "halo" and force smoothing e, a maximum red-shift 
range over which the coUisionless limit may be represented with desired accuracy. We 
also compare our results to the so-called Layzer Irvine test, showing that it provides a 
weaker, but very useful, tool to constrain the choice of numerical parameters. Finally 
we outline in some detail how these methods could be employed to test the choice of 
the numerical parameters used in a cosmological simulation. 

Key words: Galaxy: halo; Galaxy: formation; globular clusters: general; (cosmology:) 
dark matter; (cosmology:) large-scale structure of Universe; galaxies: formation 



1 INTRODUCTION 

Numerical simulations of structure formation in the uni- 
verse in cosmology use the A'^ body method in which the 
continuum density field of dark matter is represented by a 
finite number of discrete particles interacting by a smoothed 
Newtonian two body potential. It is evidently of impor- 
tance to control as much as possible for their precision and 
reliability. Specifically, beyond issues of numerical conver- 
gence, it is important to understand the limits imposed 
on the accuracy of results by the use of a finite number 
of particles to represent the theoretical continuum den- 
sity field, and the associated introduction of a smoothing 
scale (or equivalent) in the gravitational force. This lat- 
ter scale, e, clearly imposes a lower limit on the spatial 
resolution, so in order to optimize resolution the question 
is how small a value of e may be employed for a given 
number of particles and starting redshift. This question has 
been the subject of some controversy, notably concerning 
whether values of e smaller than t he initial inter-particl e 
distance may be employed (see e.g. ISplinter et al.l (|l998h : 



iKnebe et all (|2000[) 



( 2005h : |jovce et al 



Power et all (|2003| \: iHeitmann et al.l 
20081 ) : lRomeo et all (|2008t )'). 



In this article we discuss one way in which cosmological 
N body codes may be tested for their reliability which has 
not been explored previously. The idea is based on the simple 
observation that, applied to the simulation of an "isolated" 
overdensity (i.e. a finite system of size much smaller than 
that of the periodic simulation box), a cosmological simu- 
lation should be equivalent, in physical coordinates, to one 
performed in open boundary conditions without cosmologi- 
cal expansion. Indeed the only differences between the two 
should arise from possible differences in the force smoothing 
and finite size effects, both of which are variables on which 
the physical results of a cosmological simulation should not 
depend. Even without a direct comparison with simulations 
in open boundary conditions, the desired coUisionless na- 
ture of cosmological simulations can be tested for by prob- 
ing whether an isolated virialized structure, corresponding 
to a coUisionless equilibrium, remains stable in physical co- 
ordinates. 
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By "isolated" we mean that there is no other mass in 
the periodic box other than the structure considered, which 
itself evolves in a region of a size small compared to that 
of the box. The structure is therefore isolated but for the 
interaction with its "copies" included in the infinite system 
over which the force is summed. We illustrate with a set of 
numerical simulations how this required equivalence of the 
evolution in codes with and without expansion can be used 
to actually determine whether a given choice of numerical 
parameters for cosmological simulations is appropriate. We 
focus in particular on the choice of the smoothing length in 
the force, and show that the test allows one to determine a 
range of appropriate values. 

To avoid possible confusion it is probably useful to un- 
derline the distinction between stable clustering as we study 
it here, and the same term as it i s frequently discussed i n 
cosmological simul ations (see, e.g.. lEfstathiou et al.l (flQSSh : 
ISmith et al.l l)2003l ')l: it can be postulated (jPeebled [l980l ) 
that, in the strongly non-linear regime, structures evolve 
OS if they were isolated from the rest of the mass in the 
universe. If this "stable clustering hypothesis" is valid (to 
a good approximation, on average) it leads, when matched 
with linear theory, to very specific predictions for the nature 
of the correlations in the non-linear regim^H- Here, in con- 
trast, we will consider by construction the evolution only of 
a single structure, for which stable clustering must be ob- 
served if the simulation is reproducing the desired physical 
limit. 

The article is organized as follows. In the next section 
we show explicitly how cosmological codes used to simulate 
isolated overdensities in an expanding universe are related, 
in an appropriate limit, to non-expanding codes. We then 
discuss the particular limit of virialized halos, for which the 
stationarity in a non-expanding code corresponds to stable 
clustering in physical coordinates in the expanding code. In 
the following section we illustrate and investigate the test 
using a set of simulations in an Einstein de Sitter (EdS) uni- 
verse, which differ only in the smoothing length employed. 
We then also consider a set of simulations in which only the 
box size is varied. In the subsequent section we compare the 
test with the so-called Layzer Irvine test for the evolution of 
the energy in cosmological simulations. In Sect. [5] we discuss 
what can be inferred from our results about the role of dif- 
ferent possible discreteness effects in producing the observed 
deviations from stable clustering, and what can be concluded 
about dependences of these deviations on the relevant pa- 
rameters (particle number, force smoothing, box size). In the 
following section we specify in a "recipe" form how the test 
could be employed practically by simulators to test choice of 
numerical parameters in cosmological simulations. We con- 
clude with a brief discussion of possible variants on the tests 
and some more general comments. 



^ Specifically starting from power-law initial conditions it leads 
to the prediction of a "stable clustering hierarchy" characterized 
by a power-l aw correlation function of which the exponent can be 
determined l |Peeblej|l98ol) 



2 EVOLUTION OF AN ISOLATED 

OVERDENSITY IN COSMOLOGICAL N 
BODY CODES 

Dissipationless c osmolo gi cal N-body sim u lation s ( see, e.g., 
iKravtsov et"al] (|l997l) ; ISpringel et aD (|200ll ): iTevssieij 



(|2002h , 

equations 



1 ■ ■ — 

Baglal (|2005f ) for a review) solve numerically the 



+ 2H 



dt 



where 



Xi 



X,; 



(1) 



(2) 



Xi are the comoving particle coordinates of the i = 1...A*' 
particles of equal mass m, enclosed in a cubic box of side 
L, and subject to periodic boundary conditions, a{t) is the 
appropriate scale factor for the cosmology considered, and 
H{t) — d/ais the Hubble constant. The function We is a reg- 
ularisation of the divergence of the force at zero separation 
— below a characteristic scale, e, which is typically fixed in 
comoving units. For simplicity we will drop this function in 
this section as it plays no role for our considerations here. 

The superscript 'P' in the sum in ([2]) indicates that it 
runs over the infinite periodic system, i.e., the force on a 
particle is that due to the A*' — 1 other particles and all their 
copies. The sum, as written, is formally divergent, but it is 
implicitly regularized by the subtraction of the contribution 
of the mean mass density. This is physically appropriate, in 
an expanding universe , as the mean mass density sources 
the expansion (see e.g. IPeeblesI (llOSd )). 

2.1 Equations of motion a single "isolated" 
structure in physical coordinates 

Let us now consider the case illustrated schematically in 
Fig[T]where the N particles are contained within a spherical 
volume, n, of radius R, with R < L/i (where L is the side 
of the cube) . The latter condition is sufficient to ensure that 
the distance between any particle i and any other particle j 
in is less than that separating i from any copy of j in the 
infinite periodic system. The force on a particle i may then 
be written 



F, = Ff + F^ 



where 



F ' 



-GmE 



Xi 



Xi 



(3) 



(4) 



is simply the direct one over the N — 1 other particles in the 
volume SI, and 

Xi — X, — nL 



Xj 



(5) 



X, - nL13 

is the force due to all the particles in the copies, labelled 
by all vectors of non-zero integers n, and where we have 
|xi— Xj I < L/2 < \nL\ . Note that Fp is clearly finite and well 
defined, while each sum over n giving Fi^ is now formally 
divergent, but again implicitly regulated by the subtraction 
of the mean density. 
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Vi = a(t)xi, simply as 



-GmY,- 



(9) 





Figure 1. Schematic representation of the case studied in this 
paper: a single structure of A'^ particles confined in a region of 
characteristic size small compared to the side of the cubic cosmo- 
logical simulation box (dark line). 



To calculate Ff^ we observe that it corresponds simply 
to the force on a single particle displaced by (xi — Xj) off an 
infinite perfect la ttice of lattice spacing L. It is straightfor- 
ward to show fsee lGabrielli et al] (|2006l )) that, expanding in 



Taylor series in (xi — Xj) about x^ — Xj = 0, we have 
47rGpo 



2,N 



-(xi — Xj) + 0(|xi — Xj| /L 



(6) 



where po = mN/L^ is the total mean mass density (and 
thus po/N the mean mass density of the lattice of the par- 
ticle j and its copies). The leading non-zero "dipole" term 
on the right hand side in ((S)) is a repulsive term which arises 
from the subtraction of the mean mass density to regulate 
the sum: it is precisely the force arising from the mass con- 
tained in a sphere of radius |xi — Xj| of constant mass density 
—po/N. We do not write explicitly the sums for the subse- 
quent (quadrupole and higher multipole) terms, but they are 
manifestly convergent and suppressed by positive powers of 
(R/L) compared to the dipole term. 

As the sum over Xj in F^ vanishes because we have 
chosen (without loss of generality) to place the center of 
mass of the A'^ particles in f2 at the origin of our coordinates, 
retaining only the dipole term in F^ the equations of motion 
([T)) we obtain 



dxi Gm^-^ 



, 47rGpo 
H ^ — X, 



(7) 



1 particles in fi. 



where the sum is now only over the A'^ 
Assuming an EdS cosmology, for which 

d 4ttGpo 
a 3a'' ' 

these equations ([2]) may be written, in physical coordinates 



(8) 



i.e. as the equations of motion of A' isolated purely self- 
gravitating particles 0. 

Thus, when a cosmological code is used to simulate an 
isolated structure in an expanding universe, it should re- 
produce the same result, in physical coordinates, as that 
obtained for such a structure in open boundary conditions 
without expansion. This identity is valid up to 

• finite size corrections, arising from the use of a finite 
(periodic) box in cosmological simulations (and which vanish 
in the limit R/L — >■ oo). 

• eventual difi'erences due to force softening in the two 
type of codes (which we have neglected above). 

Any dependence of results on the box size or force 
smoothing is unphysical in cosmological codes. Thus these 
codes can be tested by using them to simulate isolated struc- 
tures and comparing the results obtained to those for the 
same initial conditions in open boundary conditions and 
without expansion. 



2.2 Virialization and stable clustering 

Results for the detailed evolution from arbitrary initial con- 
ditions in open boundary conditions and without expansion 
can be obtained in general only from numerical simulation. 
However, even without performing such simulations, it is 
possible to do tests of cosmological codes which are based 
on well established generic features of the evolution in open 
boundary conditions. One such feature, for a very wide class 
of simple initial conditions, is the evolution t o a virial equi- 
librium in a few dynamical times (see e.g. iHeggie fc HutI 
These equilibria are, in the limit A' — >■ oo, station- 
ary states corresponding to time independent solutions of 
the coUisionless Boltzmann equation. In a finite A' system 
they evolve away from this equilibrium, on a very long time 
scale diverging with N, due to coUisional effects. The sta- 
tionarity of these states in the coUisionless limit corresponds, 
in a cosmological code, to "stable clustering" of the virial- 
ized system. 

It follows that by studying the stability of the evolution 
obtained from appropriate initial conditions, we can test cos- 
mological codes both for effects arising from the finite size 
of the box, and force smoothing, as well as for coUisional 
effects. These are precisely the principle undesired effects 
introduced by using the N-body method to solve the (con- 
tinuum, infinite system) cosmological problem of formation 
of structure, and present even in the idealized limit that 
the numerical integration of the equations of motion is ar- 
bitrarily accurate. Thus by testing for the validity of stable 
clustering in such a regime we can test for the capacity of 



^ In the case of a cosmology with matter and a cosmological con- 
stant A, we obtain in physical coordinates an additional repulsive 
term arising from A. This can easily be incorporated in the con- 
siderations which follow, but we consider here for simplicity the 
case A = 0, i.e., the EdS cosmology. 
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where vl = a{t)-^ is the particle peculiar velocity, i.e.. 




Figure 2. Density profile (and, inset, virial ratio) at indicated 
times for our cliosen initial conditions when evolved in open 
boundary conditions (and without expansion). 



N body codes to reproduce correctly the relevant physical 
limit. It is such a test which we will focus on in what follows. 



3 NUMERICAL STUDY OF THE TEST 

To illustrate the t est we use the GADGET code 
ISpringel et al.ll2005^ . which can be used to perform both 
cosmological simulations, and simulations without expan- 
sion (in both open and periodic boundary conditions). For 
our cosmological simulations we consider, for simplicity, evo- 
lution in an Einstein de Sitter background. All our simula- 
tions here are for A'^ = 10* particles. 



3.1 Initial conditions and choice of units 

The initial conditions we study here are the following: the 
N — 10* particles are randomly distributed in a spherical 
volume of radius R = O.IL, and assigned random velocities 
sampled from a probability distribution which is uniform in 
a cube centered at zero velocity. These velocities are then 
normalized so that b — —1, where 



is the virial ratio, and K, 



^ 2Kp 

Up 

and Ur 



(10) 

are the peculiar kinetic 



energy and peculiar potential energy respectively. These are 
defined by 



(11) 



and 



Up = ^ ^diWi - Xj\) 



(12) 
(13) 



where g{r) is the exact (GADGET) two body potential. 
Thus, modulo force smoothing, Up is equal to the Newto- 
nian potential energy in physical coordinates, and we will 
therefore refer to it as the physical potential energy^ 

Our motivation for this choice is that it is a simple one 
which, although out of equilibrium, rapidly settles to a virial 
equilibrium. We note that it corresponds, in the cosmological 
simulations, to an initial density fluctuation of amplitude 



5=^ = 
po 



_3£_ 



240 



(14) 



where po is the mean (comoving) mass density of the uni- 
verse. Thus it can be thought of, roughly, as representing an 
almost virialized spherical halo at its formation time, which 
is then evolved forward in isolation from the rest of the mass 
in the universe 0. In our conclusions we will briefly discuss 
other initial conditions which it would be relevant to study 
in testing cosmological simulations. 

The results we report require only choice of units for 
length and energy: for the former we will take units defined 
by L = 1, and for the latter units in which {3GM^ /5R) — 1, 
i.e. in which the initial continuum limit potential energy is 
(minus) unity. 

We note that p4p implies that, at expansion factor a, 
starting from a = 1 ai t = to = l/^6nGpo, we have 



t — to — [d 



3/2 



See |http : //miv . mpa-garching . mpg . de/gadget / 1 



l]r,c ~ 6.6 [a^'^ - 1]tsc 
(15) 

where t^c = \/3n /32Gp is the collapse time for a cold uni- 
form sphere with mass density p. 

Our expanding simulations are evolved up to a scale 
factor a = 20, unless otherwise indicated. This means our 
study is (roughly) of halos of A = 10* particles which formed 
at a red-shift less than about twenty. Note that (|15|l implies 
that a = 20 corresponds to several hundred dynamical times 
of the halo. For N = 10* particles this is sufficiently long, 
as we will see, to see evidence large deviations from stable 
clustering in all our simulations. 



3.2 Parameters of expanding simulations 

We consider a set of simulations (of N — 10* particles) with 
the values of smoothing e shown in Table [T] In GADGET 



* Note that g{r) differs from the exact two body potential used 
in the dynamical evolution, because of the modifications associ- 
ated with the periodic boundary conditions; often the nomination 
"peculiar potential energy" is used for Up defined as in Eq. I I13I I 
but including this modification in g{r). 

^ The virial condition 6 = — 1 as given indeed corresponds, to a 
very good approximation, to the more evident condition in phys- 
ical coordinates: 2Kp = —Up implies ~ GM/R, from which 
it follows that v^/H'^r? 5, i.e., the peculiar velocity is large 
compared to the Hubble fiow velocity. 
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Name 




e/L 




e/e 


si 


3.7 


X 10- 


-5 


0.0008 


s2 


5.0 


X 10- 


-5 


0.001 


s3 


8.0 


X 10- 


'5 


0.0016 


s4 


1.0 


X 10" 


'4 


0.002 


s5 


3.7 


X 10" 


-4 


0.008 


s6 


5.0 


X 10" 


-4 


0.01 


s7 


3.7 


X 10" 


-3 


0.08 



Table 1. Names and corresponding values of e of simulations 
with N = 10* particles. L is the box size and i = LN~^/^ . 



the smoothed two body potential has a complicated func- 
tional form which is a spline interpolation between the exact 
Newtonian potential — above a separation of 2.8£ — and a 
potential of which the first derivative vanishes at zero sepa- 
ration. The value of e we quote is the value of the parameter 
with this name in the code. At this separation the smoothed 
potential is approximately 75% of it Newtonian value, while 
at e/2 it is down to approximately 50%. 

The values of the parameters controlling precision of the 
time stepping and force calculations are given in Appendix 
\K\ as well as a discussion of tests we have performed for the 
sensitivity of our results to variation of these parameters. 

Each row in Table [T] gives the name of a simulation and 
the corresponding value of e in units of the box size L. Also 
given is the ratio of e to f = LN~^^^ . The latter corresponds 
to the initial grid spacing of a cosmological simulation with 
the same mean (comoving) matter density. The initial mean 
nearest neighbour separation A, on the other hand, is given 
byl!l 



A = 0.55 X 



3iV 



1/3 



0.9^ 



R 



(16) 



i.e. e/A « Qe/l since R = O.li. 

We thus consider in all our simulations, as in many large 
cosmological simulations, a smoothing which is fixed in co- 
moving units. As we will discuss a little more in our con- 
clusions, our test can of course be applied with any other 
smoothing prescription, e.g., fixed smoothing in physical co- 
ordinates or adaptative smoothing. The motivation for the 
range of e we have chosen to study, which extends to values 
significantly smaller to those typically used in large cosmo- 
logical simulations, is the following: 

• An upper cut-off on e is imposed by the fact that this 
scale must be sufficiently small so that gravitational mean 
field forces may be well approximated. This requires simply 
that 



e<s:Rs 



(17) 



where Ra is the characteristic scale on which the mass den- 
sity in the structure varies (in the continuum limit). The 
value of e in s7 corresponds to e ~ 0.04 R. Given that e is 
fixed in comoving coordinates, while stable clustering will 
lead to a structure with Rs oc 1/a, this simulation (run up 



^ This expression is derived from the value for an infinite Poisson 
distribution (for which the exact neare st neighbour distribution 
is known, see e.g. iGabrielli et al. 



to a = 20) will clearly be expected to manifest the effects of 
the violation of the bound (|17p . 

We note that even the largest value of e we consider corre- 
sponds to a value significantly smaller than I. Indeed such a 
choice is unavoidable if one wishes to resolve the non-linear 
evolution of structures with a modest number of particles in 
a cosmological simulation. As mentioned in our introduction, 
the use of e < ^ has been the source of discussion and con- 
troversy in the literature, notably as at early times it leads 
inevitably to effects whic h should be absent i n the desired 
mean-fi eld limit (see e.g. Splinter et al.l (| 19981 ) ; I Joyce et al.l 
((200ab : iRomeo et al.1 (|2008l ')'). Our present test, which con- 
siders only the strongly non-linear regime subsequent to col- 
lapse and virialization, clearly cannot give us any informa- 
tion or constraint on the accuracy with which coUisionless 
behaviour is reproduced in the early time evolution. We will 
return briefly to these issues in our conclusions. 

• A lower cut-off on the other hand is dictated in principle 
only by numerical limitations: without such a cut-off hard 
two body collisions will occur, with (arbitarily) large acceler- 
ations requiring integration wit h correspondingly sm all time 
steps (for a discussion, see e.g. iKnebe et al.l (|200d )). Given 
that a two body collision is soft if (Gm/sv1) <^ 1, where 
Vr is the initial relative velocity and s the impact factor, a 
naive estimate of the condition on softening to suppress hard 
collisions is e 3> {Gm/v^), where v is the typical speed of a 
particle in the system. For a virialized system of A'^ particles 
of size Rs we have Nv^ ~ GmN^ / Rs- Thus we estimate 



e > Rs/N . 



(18) 



The value of e in si corresponds approximately to this esti- 
mated lower bound at the beginning of the simulation. Note 
that, in contrast to the upper bound Eq. (|17[) . an evolution 
corresponding to stable clustering (with Rs ~ 1/a) will im- 
prove progressively the satisfaction of the condition Eq. (|18p . 



Given that the goal of N body simulations is to repro- 
duce the coUisionless limit, in which hard two body collisions 
should clearly play no role, the imposition of the lower bound 
(|18|) is clearly justified. Indeed we note that simple estimates 
of the minimal e required to reproduce the coUisionless limit 
suggest that much larger values of e may be necessary. If one 
assumes, for example, that e should be large enough so that 
the force from a single particle is always subdominant with 
respect to the mean field force, one obtains 

e>i?»^f"'/^ (19) 

An even more restrictive condition, that might possibly be 
relevant, is that e should be at least of order the average 
interparticle distance, i.e., 

-1/3 



e > RsN- 



(20) 



As we will discuss further below, although we do not do 
so here the test we develop could in principle be used to 
determine which (if any) of these scalings is the right one 
for the minimal e. 



3.3 Evolution in open boundary conditions 

As discussed above, the test we explore here, for stability 
in physical coordinates of a virialized structure, does not 
necessarily require direct comparison with the same initial 
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conditions evolved in open boundary conditions. Such com- 
parison constitutes, however, a more stringent test, and, as 
we will see, can be used to derive more precise quantitative 
conclusions about the choice of numerical parameters. 

We have therefore evolved the initial conditions de- 
scribed above in open boundary conditions without expan- 
sionQ The values of the numerical parameters we have used, 
and tests we have performed for the stability of these results, 
are given in Appendix [X] Shown in Fig. [2] is the initial and 
evolved density profile obtained at the indicated times in 
a simulation in open boundary conditions and without ex- 
pansion. In inset is shown the evolution of the virial ratio 
6 (defined as Eq. pO|) . taking a{t) = I). After a few dy- 
namical times the system has settled down "gently" to a 
macroscopically stable virialized configuration, with small 
residual fluctuations about it. Below we will use a smooth 
fit to the profile at f2rsc as a template for the profile of 
the collisionless equilibrium established in open boundary 
conditions from these initial conditions. We will then com- 
pare this template with the profile of the virial equilibrium 
obtained in the cosmological code. In absence of coUisional 
effects, finite size effects or other numerical effects, the two 
should coincide. Note that one could also perform a test in 
which the open system is evolved to the considerably longer 
times on which coUisional effects play a role, and compare 
this with the evolution in the expanding case. This, however, 
is not a test we explore here as our focus is on testing the 
collisionless nature of cosmological simulations. 

3.4 Evolution in cosmological code 

3.4-1 Virial ratio 

In Fig. O are shown for the indicated simulations, the evolu- 
tion, as a function of redshift, of the virial ratio, as defined 
by Eqs. dni) and 

For all the simulations, except s7, the virial ratio evolves 
qualitatively as would be expected if the system evolves as 
in open boundary conditions: they show low amplitude co- 
herent oscillations which decay gradually indicating viriali- 
sation (corresponding to b = —1). Further we have checked 
that there is good quantitative coherence between the ampli- 
tude and time scale of these oscillations, and those found in 
open boundary conditions (inset of Fig. [5]): using Eq. (|f 5p 
we see that the temporal range of the latter corresponds 
just to evolution up to a ~ 2 in the expanding case. We will 
consider below the exact degree of agreement between the 
density profiles obtained in the expanding simulations and 
the non-expanding case. 

The fact that s7 behaves so differently — deviating 
clearly from a behavior like that in the other cases at a ~ 5 
— can be attributed, as anticipated above, to the violation 
of the constraint (|17[l: we have e/Rs ~ 0.04 initially, which 
means that, if the structure does indeed remain fixed in 
physical coordinates, at a ~ 5 the effective mean-field force 
due to all particles is very different to its Newtonian value. 
For s6, on the other hand, with a smoothing about seven 
times smaller than in s7, we expect to have e/Ra « O.I at 

^ See, e.g., IWorrakitpoonponI l l201ll) : ISvlos Labinil | |2012|) for a 

more detailed study of this case. 




■ 1 10 

a 

Figure 3. Evolution of the virial ratio b (as defined in the text) 
up to scale factor a = 20, for the different simulations in Table [T] 

a = 20, still small enough so that the Newtonian value of the 
mean field potential is well approximated by the smoothed 
potential. In this latter case, however, even from a ~ 10 the 
virial ratio already appears in Figs [2] and [3] to show a slight 
tendency to increase at the larger values of a. 

3.4-2 Potential energy 

In Fig. 2] are shown the evolution of the physical potential 
energy Up. If the structure, once virialized, remains stable 
in physical coordinates we should have Up = constant. Fur- 
ther, of course, if the smoothing plays no role, this constant 
value should be the same in all the simulations. We observe 
that only the behavior observed for the simulations s5 and 
s6 appears to be consistent with stable clustering of the viri- 
alized structure, and even in these two cases a slight devi- 
ation is apparent at the end of simulated red-shift range, 
from a ~ 18. These plots thus indicate that at most in the 
corresponding narrow range of e can the behavior required 
in the desired continuum limit be reproduced by the TV body 
method. 

3.4.3 Density profiles 

Let us now examine whether this conclusion is borne out 
by further analysis of the evolved configurations. Shown in 
Fig. [5] are the measured density profiles at the indicated 
scale factors in each of the simulations. The results confirm 
strongly what can be anticipated from the analysis of the po- 
tential energy: the profiles agree increasingly poorly in time, 
with s7 and si clearly giving profiles completely different to 
those obtained in the other cases. s4, on the other hand, 
shows a much smaller discrepancy with the remaining two, 
s5 and s6. These latter two simulations agree very well with 
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Figure 5. Density profiles in comoving coordinates as measured in the indicated simulations from Table [T] at the different indicated 
scale factors. 



one another, except at the very last plot where a slight dif- 
ference may be observed. We note that, compared to these 
cases, the simulations with a smaller e have a denser more 
concentrated core, while for s7, which has a larger value of 
e, the opposite behavior is observed. This is very consistent 
with our comments above about this latter simulation: e is 
so large that mean field forces are very reduced compared 
to the exact Newtonian mean field, leading to a much less 
condensed structure. 



3.4-4 Rescaled density profiles 

It can be seen qualitatively from the previous figures (which 
are plotted in comoving coordinates) that the comoving size 
of the structures does indeed decrease, with a correspond- 
ing increase in their density. To see whether the behavior is 
quantitatively in agreement with that associated with stable 
clustering, we show in Fig. [6] the evolution of the profile in 
physical coordinates for the simulations indicated, i.e., we 
plot in each case n(r) = n{x)/a?' as a function of r = ax. 
In this representation stable clustering corresponds to an in- 
variant profile. We also show in these plots the template for 
the profile of the coUisionless equilibrium obtained from a 
simulation of the open non-expanding described in 
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1 10 

a 

Figure 4. Evolution of the potential energy Up of the structure 
(in physical coordinates), for the indicated simulations in Table[Tl 



Sect. 13.31 The insets in the plots show the results in each 
expanding simulation normalized to this profile. 

The results confirm what has been anticipated above 
from the examination of the behaviour of Up and the com- 
parison of density profiles, but also give additional con- 
straints: s5 and s6 reproduce stable clustering considerably 
better than any of the other simulations, but s5 also clearly 
does better than s6, which shows a deviation between a — 15 
and a = 20. Thus, of the full range of e considered, s5 is 
closest to optimal, while all others lead to quite measurable 
deviations from the continuum limit behaviour. We see very 
clearly in these plots again the marked qualitative differ- 
ence between the cases of a smoothing which is too large 
and one which is too small. In the former case the structure 
obtained is very much less dense and more extended than 
it should be, while in the latter case it is very much more 
compact, with a density profile which steepens towards the 
centre. These differences, as we will discuss further in the 
final section below, clearly correspond to the very different 
effects at play in the two cases. 

We note that the conclusions in the previous paragraph 
can be drawn even without the direct comparison with the 
non-expanding density profile: in other words, when stabil- 
ity is observed in a given expanding universe simulation, the 
(stable) profile obtained is always consistent with the correct 
one. The comparison with the non-expanding profile, shown 
in the inset of each figure gives a more quantitative measure 
of the deviation of the corresponding expanding universe 
simulation from the correct behavior. Note that these insets 
have been cut in all cases at r = 0.1, since in all cases there 
are very significant deviations beyond this radius, which is 
reflected also clearly in clear deviations from stable cluster- 
ing. Indeed we note that in all cases the very outer part of 
the profile extends very significantly further than in open 
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Table 2. Parameter of simulations with N = 10* particles for 
various different initial system size R (in units of the side of the 
periodic box). The ratio e/R is the same as that in the simulation 
s5 in Table [T] The parameter Oond is defined in the text. 

boundary conditions. In most cases this corresponds to a 
very small fraction of the total mass, except in the case of 
s7 (with the largest smoothing) for which the characteristic 
size of the whole structure is, as we have noted, very much 
larger than it should be. 

3.5 Test for box size dependence 

The simulations in the previous section are of fixed box size. 
The fact that they evolve differently is, by construction, due 
only to the different force smoothing. However when we com- 
pare the profiles obtained to the template in open bound- 
ary conditions (determined at very early times and taken to 
be representative of coUisionless evolution), differences may 
also arise because of the periodic boundary conditions. As 
we have explained the differences with respect to the open 
system should be suppressed by powers of R/L and thus van- 
ish as the size of the system becomes small compared to the 
box size. By varying the box size we can test both whether 
finite size effects may be observed, and whether they are, as 
we have assumed, small for R/L = 0.1 

To do so we have run the set of simulations in Table [J] 
The simulation s5R0.1 is identical to the simulation s5 con- 
sidered above, except that it has been run for a longer time, 
up to a = 33 (rather than a = 20). The other four sim- 
ulations are for the same number of particles, A'^ = 10'', 
and differ only in their initial size R. Because the overden- 
sity 5 represented by these initial conditions depends on R 
[cf.Eq. (|14|) ]. the relation Eq. (|15[l between the physical time 
elapsed and the scale factor a is modified. More precisely, in 
units of the characteristic time of the isolated structure Tsc, 
a given elapsed time t — to corresponds to a fixed value of 
^-3/2 (-^3/2 _ j^-j rpj^g g^^jg factor a = Oend in Table [2] up to 

which the corresponding simulation has been run, has been 
chosen so that R~^^'^ (^^cnd ~ 1) equal in all cases, i.e., all 
simulations are run up to the same time in units of r^c. 

The choice of e given in Table [2] have been made by 
scaling the value in s5 in proportion to R (and the mean 
interparticle distance). As e is fixed in comoving coordinates, 
it evolves in physical coordinates in proportion to a(t), and 
therefore as a function of {t — to)/Tsc in a manner which 
depends on the box size. Thus the dependence on box size 
in the evolution of these systems can arise not just from 
contributions to the gravitational force due to the periodic 
copies, but also through possible differences in the dynamics 
due to the differences in force smoothing in each case. Apart 
from such effects their evolution should be identical when 
analysed in physical units (of length and time). 
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Figure 6. Density profiles in physical coordinates at the different indicated values of a. Each figure corresponds to the single indicated 
simulation. Also shown in each case is a smooth fit to the profile (labelled OBC) obtained in an open non-expanding simulation at 
t = 12tsc- In the inset plots of the profiles normalized by this latter profile arc given. 



In Fig. [7| is shown the evolution of the potential en- 
ergy Up, normalized to the modulus of its initial value, for 
the five simulations, as a function of the appropriately nor- 
malised time determined by Eq. (|15p . and the evolution of 
the virial ratio in the inset. In Fig. [8] is shown the density 
profile in each of the five simulations, at the time in each case 
corresponding to a = 10 in the simulation s5 (and s5R0.1). 

Very strong dependence on the box size is manifest for 
the two largest systems (of which the initial diameter is equal 
to, respectively, 0.4 and 0.6 of the box size): both the energy 
of the virial equilibrium attained, and the density profiles, 
are very significantly different to those in the other cases. 



The smaller systems, on the other hand, show a clear con- 
vergence for these quantities (and which we have seen for 
the case s5R0.1 agree well with those obtained in the open 
case). The differences at late times in the deviations towards 
higher values of Up, associated (as can be seen in the inset) 
with deviations from satisfaction of the virial condition, may 
clearly be attributed to the differences in the force smooth- 
ing noted above: a larger box size corresponds, at a given 
{t~to)/Tsc, to a larger scale factor a, and therefore to a larger 
smoothing with respect to the size of the system. Indeed in 
the simulation s7 we saw (Fig. [4] ) that significant deviation 
in the behaviour of Up already at a ~ 2 — 3. The simulations 
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7. Evolution of the potential energy Up for the simula- 
Table [2] with N = 10"*, normalized to initial absolute 



Figure 8. Density profiles for the simulations in Tabic [2] with 
N = 10'*, at the normalized time corresponding to a = 10 in 
S5R0.1 



here have e/Rs initially exactly ten times smaller than in s7, 
and so would be expected to show deviations in the range 
o ~ 20 — 30 as is observed. 

Thus for the quantities we have focussed on above the 
effects due to the finite box size in the system of the ini- 
tial size we have considered appear indeed to be very small. 
We note, however, that considerable box size dependence is 
manifest even for the smaller boxes in other quantities. The 
amplitude of the oscillations during the initial relaxation to 
virial equilibrium are very significantly larger than in the 
smaller systems, for which the amplitude appears to con- 
verge approximately. This suggests that the corresponding 
modes of oscillation of the structure about the virial equi- 
librium are enhanced very significantly by the gravitational 
coupling to the periodic copies. 



4 COMPARISON WITH THE LAYZER-IRVINE 
TEST 

One possible test of the accuracy of a cosmological code is 
the so-called Layzer-Irvine (LI) test, derived from the equa- 
tion of the same name which describes the variation of tot al 
energy in an expanding universe (see, e.g.. |PeebleiJ (|l980l )l. 
Unlike energy conservation in non-expanding simulations, it 
is, however, a test which is rarely employed in practice by 
cosmological simulators as a control on their code, because 
it is not evident how to quantify the violation of the LI equa- 
tion which may be toleratecjf]. Given that the test discussed 



* The GADGET2 user guide, for example, states that "the cos- 
mic energy integration is a differential equation, so a test of con- 
servation of energy in an expanding cosmos is less straightforward 



in this paper is an independent one on the correctness of 
cosmological codes in a specific regime, it can in principle 
be used to "calibrate" the LI test in this context. Likewise 
it is interesting to see whether it may be useful to employ 
both tests together. 

In terms of the quantities defined above , the Layzer 
Irvine equation may be written (|Peebleslll98"ol ) 



^ [a(AV + Up)] 
We thus define the quantity 



a{Kp + Up) + /; Kpda 



(21) 



(22) 



Kpil) + Upil) 

which should be equal to unity 0. 

In Fig. [9] is shown the evolution of A{a) for the sim- 
ulations si to s7. We observe immediately that the two 
simulations in which A{a) remains close to unity are pre- 
cisely those, s5 and s6, which have been singled out by our 
test above as reproducing best the required behaviour. On 
the contrary, all the other simulations which showed much 
greater deviation from stable clustering also show larger de- 
viations from unity of yl(a). More precisely, in all cases where 
deviation of A{a) from unity by more than a few percent is 
observed, the stable clustering test showed the results for 



that one may think, to the point that it is nearly useless for testing 
code accuracy." 

^ The LI test remains applicable in this form for the infinite peri- 
odic system, if Up is calculated with the corresponding two body 
potential. As this makes no significant difference to the results we 
give below, we continue to use Up as considered elsewhere in the 
paper (i.e. calculated without this modification to the two body 
potential) . 
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the clustering in the system were completely incorrect and 
unphysical. 

The correlation between the information we deduced 
from Fig. U and what we observe in Fig. [^is very evident, 
and the reason for it very simple: the requirement that Up 
be constant is, when the virial condition b = —1 is satis- 
fied, equivalent to the condition that both Kp, and therefore 
also the sum Up + Kp corresponding to the energy in phys- 
ical coordinates, is constant. In this case, as can be seen 
from Eqs. H21I) and H22[) . the LI test is also satisfied. Thus 
it is clear that the very strong deviations from stable clus- 
tering observed for the simulations si to s4 stem from a 
poor integration of the equations of motion, with very sig- 
nificant violation of energy conservation. As discussed in 
Sect. 13.21 such numerical difficulties are expected to arise 
due to the precision requirements of integrating accurately 
the hard two body collisions present when e becomes very 
small. We note that our results are very consistent with 
those of iKnebe et"al] (|200(]| ') who have studied in detail the 
effects of such poorly integrated hard collisions: they lead 
to an artificial injection of energy into the system, increas- 
ing its size as some particles are sent into spurious higher 
energy orbits. We observe here indeed quite distinctly these 
effects both in the behaviour of the energy (which increases) 
and the profiles which stretch out further than they should 
(compared to stable clustering). 

For the very small e considerably greater precision than 
that employed (see Appendix|X| would be required to attain 
numerical convergence. We note that our results do suggest 
that, at given numerical precision, a lower bound on e may 
be expressed in a simple form like Eqs. (|18|I - (|2U| I. i.e., that 
the minimal e required scales linearly the size of the struc- 
ture R3 : in Fig.|4]each of the simulations si to s3, which have 
the same maximum time step, show an approximate plateau 
in Up, starting from a scale factor a which increases roughly 
in proportion to 1/e. Given that we observe in these simula- 
tions that Rs ~ 1/a, this behaviour of Up, indicating energy 
conservation, therefore sets in approximately at a fixed value 
of e/Rs, in line with bounds of the form of Eqs. (fT8)) - (f20)) . 
A study of simulations with different A'^ would be required 
to establish which (if any) of the proposed scalings is the 
correct one 

We underline that the LI test for an expanding simu- 
lation is, just as the test for the constancy of Up, a weaker 
test than the test for the stability of clustering: the latter 
tests for the coUisionless nature of the evolution, which is a 
different (and stronger) requirement than energy conserva- 
tion. In practice, however, the breakdown of energy con- 
servation is often due to the difficulty of integrating nu- 
merically with sufficient accuracy the coUisional dynamics 
(specifically, hard two body collisions), and therefore the 
breakdown of the coUisionless approximation is associated 
with the violation of energy conservation. Such an associ- 
ation can always be "undone" , in principle, by increasing 
sufBciently the accuracy of the numerical integration. In 
practice, however, it is very difficult to disentangle the two 
effects, and indeed studies up to now of two body c o Uision - 
ality in cosmological simulations (e.g. iKnebe et al.l (|2000l ): 

1° iKnebe et all | |2000D propose bounds based on the same simple 
argument given above for Eq. III8I I. 
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Figure 9. Layzcr Irvine test for the simulations in Table [T] See 
text for the definition of A{a). 

IPower et all (|2003h : iBinnev fc Kneb^ (|2002l 'l'l have not done 
so. 

In conclusion we find that the LI test is a very use- 
ful and relevant one in the context of the present test of 
cosmological simulations. Quite simply deviations of the di- 
mensionless parameter A defined above by more than a few 
percent appear always to be indicative of a grossly incorrect 
evolution. The crucial difference with respect to its use for a 
full cosmological simulation, which, as mentioned, has been 
found to be problematic, arises from the difference in initial 
condition: for the very cold and almost perfectly uniform ini- 
tial conditions of cosmological simulations the denominator 
in A approaches zero, which makes it difficult to calibrate 
the test. 



5 NATURE OF DISCRETENESS EFFECTS 
AND PARAMETRIC SCALINGS 

Our results above establish that the tests considered can 
clearly detect and measure discreteness effects in A''-body 
simulations, i.e., deviations of the results of such simulations 
from the desired continuum limit. We discuss now briefly 
what conclusions may be drawn about the nature of these 
discreteness effects. More specifically we discuss what we can 
conclude about the parametric dependences of these effects. 

Discreteness effects here can be divided into two cate- 
gories as follows: 

• "numerical discreteness effects" arising arise from limi- 
tations on precision in the integration of the A^-body system 
with a given smoothed two-body potential, and 

• "physical discreteness effects" arising from the use of 
a finite particle density, finite force smoothing and a finite 
periodic box. 
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The former are related to the choice of the numerical 
parameters controlling accuracy of the force calculation and 
time stepping in the code, while the latter are related to the 
number of particles -N, the size of the force smoothing e and 
the size of the box L[^. The two kinds of effects are in prac- 
tice interrelated since the numerical precision required will 
depend typically on the values of the "physical" discreteness 
parameters. However, in order to understand the effects at 
play in cosmological simulation, it is useful to separate them 
in this way. One can then consider, on the one hand, the is- 
sue of numerical convergence at fixed values of A'', e and L, 
and, on the other hand, the scalings with A'', e and L, of the 
deviations from the desired physical behaviour, assuming 
"perfect" numerical convergence. It is the latter we consider 
here. Our results in Sect. 13.51 showed up clearly the pres- 
ence of effects related to the box size L, which, in line with 
expectations, decrease strongly as L increases compared to 
the size of the simulated structure. We do not pursue fur- 
ther tests here to establish exactly the associated scalings, 
but instead focus on the other two parameters. 

Our results show clearly, as anticipated, that for larger 
values of the smoothing, deviation from the coUisionless self- 
gravitating limit arises predominantly from the associated 
loss of spatial resolution. This was most easily "diagnosed" 
by the behaviour of the virial ratio, which deviated clearly 
away from b = —1 towards less negative values. As dis- 
cussed in Sect. 13.51 the behaviours observed are very con- 
sistent with the simple bound Eq. (|17p . with significant de- 
viations becoming easily visible (in potential energy and 
profiles) roughly when (e/Rs) ~ 0.1. Using the fact that 
Rs oc 1/a, we deduce that the scale factor at which we see 
the effects of the finite resolution start to significantly mod- 
ify the structure is 



10 



(23) 



where R is the size of the structure at a = 1 (and the result 
is valid for the case we consider where e remains fixed in 
comoving coordinates). 

The very large deviations from stable clustering we have 
observed in our simulations with very small e are, as we have 
discussed, apparently due to poor integration of hard two 
body collisions. This problem is clearly diagnosed very well 
using analysis of the Up, or even more clearly using the LI 
test. As mentioned, such effects h ave been dia g nosed and 
discussed in some detail notably in iKnebe et aP lj2000h . 

In principle, as we have discussed, two body coUi- 
sionality, when integrated accurately, can also contribute 
to the deviations from stable clustering we observe. In- 
deed in our simulations s5 and s6, which satisfy quite pre- 
cisely the LI test, we see clear deviations which are sim- 
ilar qualitatively to those observed in the case of poorly 
integrated collisions: tails in the density profiles which be- 
come more extended in time and a hint of steepening of 
the inner density profile. It is straightforward in our case 



We do not consider here the starting red-shift Zi for a simula- 
tion which is a parameter introduced in the A'^-body d iscretisation 
and on which discret e ness e ffects may depend (see e.g. ljovce et al.l 
ll2008l '); [Knebe et al.l ||2009| ')). As we study here the evolution only 
of non-linear structures from the time they form, we cannot con- 
strain Zi which can affect the evolution prior to this time. 



to estimate the time scale for such two body effects us- 
ing the well known results for the case of an open virial- 
ized system in a no n-expanding space. I n this case n umeri- 
cal stu d ies (see e.g. iFarouki fc Salpeteij (Il982l. Il994l ') : lThei^ 



J 19981): iTheis fc SpurzemI (| 19991 ): iDiemand et al. 



iGabrieUi et all (j2010l )) have shown that the time scales for 



evolution of coUisionless equilibria is very consistent with the 
tho se estimated analytica lly (originally by Chandrasekhar, 
see I Chandrasekhar! (j 19431 )) for two body collisions, given by 



T2body ~ kNTc 



(24) 



where Tc is a characteristic crossing time for the structure, 
and ft is a numerical factor incorporating the "Coulomb log- 
arithm" Modulo box size effects, which we have seen are 
small, the only difference between these open systems and 
the one we are studying should arise from the difference in 
the smoothing, which in our cosmological code simulations 
is fixed in comoving coordinates and therefore varies in time 
(increases) in physical coordinates. The two body relaxation 
time is, however, only logarithmically sensitive to the lower 
cut-off in the two body interaction, and temporal variation 
of the smoothing e will lead therefore to modification of the 
numerical prefactor k in the calculation, with at most very 
weak dependence on A'^. Using Eq. 1)15^ with Eq. (|24p we 
thus estimate that two body relaxation will start to cause 
deviations in our numerically well converged test simulations 
(i.e. s5 and s6 above) when 



0.2body 



V5 



-1 2/3 



N 



(25) 



where we have taken Tc ~ Tsc- 

Values of k of order those measured in open simulations 
(with smoothing fixed in physical coordinates), typically a 
little smaller than unity, thus give an estimate for a2body 
quite consistent with the hypothesis that two body coUi- 
sionality should account for the observed deviations from 
stable clustering in the simulations in which the energy is 
conserved well. We note that these re sults also appear t o 
be very consisten t with studies such as IPower et all (|2003h : 
iBinnev fc Knebd (|2002h which have detected such effects in 
cosmological simulations. Clearly however a full study of 
the A'^-dependence of the results of our test, which we do 
not undertake here, would be required to establish whether 
the scaling Eq. (|25|l is indeed observed. It would be instruc- 
tive to couple such a study also to one including analysis 
with other indicators of coUisionality (e.g. m easures of diffu- 



sion in velocity space like those employed in IDiemand et al 



(|2004|), or using two mass species as in iBinnev fc Knebg 
(2002)). 

An important practical question in numerical simula- 
tion is whether there is an optimal value of force smooth- 
ing in cosmological N body simulations. The question may 
be posed either with respect to some set of numerical con- 
straints (quantified e.g. in terms of bounds on the numerical 
parameters controlling accuracy of integration), or in ab- 
straction from such limitations. The relations Eq. (|25p and 



This logarithmic factor is simply log(iJ/£) when e is larger 
than the minimal impact factor required for the validity of the 
soft collision approximation. This calculation is for the case of a 
time independent smoothing. 
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Eq. H23p can be clearly combined, in principle, to provide a 
prescription for an optimal value of the latter kind: although 
we have not determined explicitly, it is clear that K{e) is a 
monotonically (a priori logarithmically) increasing function 
of e, while the coefficient in Eq. (|23p is monotonically de- 
creasing. Thus by varying e at fixed an optimal value can 
be found which maximizes the scale factor at which discrete- 
ness effects modify the evolution. A much more extensive 
study, notably of A'' dependence, would clearly be neces- 
sary to determine the scaling with of such an optimum. 
It would be interesting to determin e, in par t icular , whether 
it tur ns out to be that derived by iMerrittI (| 19961 ): iDehnenI 
from simple considerations of the optimization of the 
representation of the continuum (mean-field, Newtonian) 
forces. It would evidently be interesting also to compare with 
optimization criteria derived in var ious numerical conver - 
gence studie s in th e literature (e.g. iDiemand et al.l (|2004 ): 
IPower et al.] (|2003h ). 



6 PRACTICAL GUIDE TO USE OF THE TEST 

The numerical study we have reported establishes that the 
test considered can provide non-trivial constraints on the 
accuracy of cosmological simulations. To facilitate the ex- 
ploitation of the test in practice by cosmological simulators 
we now summarize in a recipe form how it could be im- 
plemented. The results of such tests will, as we discuss, 
inevitably depend strongly on the details of the particu- 
lar simulations (cosmological model, time and length scales 
probed, quantities of interest and desired level of precision 
etc). We therefore do not attempt to summarize the infor- 
mation which can be obtained from the test in some simple 
set of rules. Instead we propose to simulators an instrument 
they can use to assess themselves the reliability of the results 
of their codes. 

Let us suppose a simulator intends to run a (dissipa- 
tionless) cosmological simulation, which has A*'o particles of 
mass m in a box of side Lq. An initial perturbation spectrum 
P{k) is given at the starting red-shift Zi, and a cosmological 
model specifying the scale- factor (e.g. A CDM). Based on 
criteria at his disposal, he chooses a set of numerical param- 
eters (time step and force accuracy criteria, smoothing e). 
The test we have considered for the accuracy of simulation 
of halos containing A'^ particles (and of mass M = mN) can 
be generalized as follows: 

1. Consider a simulation identical for all numerical param- 
eters to the full simulation, except for the box size, which is 
rescaled so that 

L = Lo(-) (26) 

This condition means that the mean matter density of the 
universe is equal to that in the full simulation, and therefore, 
in particular, the ratio e/£ has the same value as in the full 
simulation. 

2. Distribute the A'' particles in a spherical region of radius 
R = L(3po/47rpi,)^'''' where pv is the initial density of the 
halo. For py ~ 200/9o this corresponds to i? ~ 0.1 (as in our 
simulations sl-s7 above). 

3. Assign velocities with some simple velocity distribution 



(e.g. uniform or gaussian) and normalize them so that the 
initial virial ratio is unity. 

4. Run the cosmological code (with the chosen numerical 
parameters) starting from Zy{M), the (estimated) maximum 
redshift for virialization of a halo of mass M in the model 
studied. 

5. Check that the potential energy Up and virial ratio b 
display the expected physical behaviour (low amplitude de- 
caying oscillations about a constant value). Systematic de- 
viation of the virial ratio from 6 = — 1 towards less negative 
values indicates that the upper resolution on smoothing is 
violated. Plot also the parameter A{a) for the LI test. If 
deviations from unity are at the level of more than a few 
percent, the parameter choices are inappropriate to simu- 
late over the corresponding time-scale. 

6. To obtain more precise limits on the precision of halo 
profiles (in the cases where the previous tests are reasonably 
well satisfied) apply first the test for stability of clustering. 
Further quantitative limits can be obtained by comparison 
with the equilibrium profile obtained from a high-resolution 
simulation of the same initial conditions in open boundary 
conditions. 

Our simulations sl-s7 reported above test the case 
N = 10^ for a range of values of e (and other numerical 
parameters) up to red-shift of about twenty. Following the 
prescription above we would exclude all but s5 and s6 using 
the analysis of the energies and LI test beyond a red-shift of 
a few (cf. Fig. |4] and Fig. |9} . The conclusion of our further 
analysis (step 6) above was that significant departures from 
stable clustering (and from the equilibrium template) were 
observable, increasing as a function of redshift. Beyond a 
redshift of order ten, we could conclude, for example, that 
50% precision on profiles of halos over two decades in scale is 
not attainable with A'' = 10**. To determine how many par- 
ticles (and what numerical parameters) would lead to such 
a level of precision could be determined by performing the 
test for larger TV. 

Several variants of the test as described above could 
provide further or more precise constraints: 

• Different initial conditions could be used in steps 2 and 
3. In particular rather than a structure close to, but not at, 
virial equilibrium, one could take instead a structure which 
is already at equilibrium. Such an initial condition could be 
prepared numeric ally, or set up directly u s ing the Eddington 
formula (see, e.g.. iMuldrew et al.l (|201ll ): iKazantzidis et all 
(|2004h ). It would be interesting notably to study equilib- 
rium halos with profiles like those typically measured in 
simulations (e.g. NFW halos). Alternatively initial condi- 
tions could be obtained directly from the full cosmological 
simulation, by extracting typical halos at about the time 
they virialize. Techniques to do this have been developed in 
the context of convergence studies in which individual halos 
are identified and resimulated at higher resolution(see, e.g. 
I Power et al.l (|2003h and references therein). 

• To separate out possible effects arising from the box 
size, the test could be repeated in boxes of different sizes, 
and specifically in a box of side Lq. These are important 
in particular when comparison with a simulation in open 
boundary conditions is made. As the change in the box size 
leads to a change in the mean density (by a factor of N/No 
for the case of a box of size Lo), this must be accounted for 
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in applying the test. As discussed in Sect. 13.51 in an EdS 
cosmology this can done, modulo effects due to the force 
smoothing, by a simple rescaling of the time. If the role of 
smoothing is the focus of the test, and/or the cosmology 
employed is non-scale free (e.g. ACDM,) the code would 
need to be modified "by hand" to impose the evolution of the 
scale factor corresponding to the cosmological simulation. 

• Another variant of the test would be, using further the 
framework of the spherical collapse model for halo forma- 
tion, to start from the time of turn-around. At this time the 
spherical overdensity is already sufficiently large (~ 5 — 6) 
that one should be able to consider the system to a good 
approximation as isolated. The initial condition would thus 
be taken as a uniform or quasi-uniform sphere with physical 
velocities equal to zero. Given that the scale R in step 1 
would be comparable to the box size, it would be appropri- 
ate also to study simulations in a larger box as discussed in 
the previous paragraph. 



7 CONCLUSIONS AND DISCUSSION 

The central point of this paper is the introduction of a sim- 
ple test on the reliability of cosmological A^-body simulations 
in the non-linear regime. We have shown with a very simple 
implementation of the test that it can constrain strongly the 
choice of the unphysical parameters introduced by the N- 
body method, given desired precision criteria on the prop- 
erties of virialized structures. More specifically we have illus- 
trated that the test can determine a window for an appro- 
priate force smoothing e for simulations with a given time 
stepping accuracy, and any given A in a virialized struc- 
ture. At larger e the loss of resolution was clearly seen to be 
the dominant effect, while at small e the effects of two body 
coUisionality become the primary cause of deviation. In con- 
trast to most other methods explored previously m the liter- 
ature these effects are detected and quantified by comparing 
the evolution of the test simulation with an exact behaviour, 
rather than by a comparison between different codes. 

We have not attempted here, and indeed it is not our 
goal, to try to use the test to derive some set of simple 
short-hand rules that could be used to, say, choose e in a 
given cosmological simulation (or compared with other pre- 
scriptions which have been given in the literature). Rather 
it is the simulator who should apply the test to derive what 
his constraints are, with respect to his own precision re- 
quirements. Nevertheless it is interesting to comment a little 
more on what our specific test, of simulations with fixed co- 
moving softening, indicates. Smoothing in simulations with 
fixed comoving softening are characterised by the sole pa- 
rameter e/i. We can compare the smoothings we have con- 
sidered (Table [l| to typical values employed in large volume 
cosmological simulations: the smoothing in s5, s6 and s7 
corr espond, respectively, to e/i « 1/125, 1/100 . 1/12, while 
e.g..lSpringel et al.1 (|2005h uses e/i ^ 1/50, and lSmith et al.1 
(l20oir e/i « 1/16. From the results discussed above it is 
clear that application of the test for these specific values can 
give very strong indications of the reliability and /or preci- 
sion of the properties of halos in such simulations. Indeed 
we note that, using Eq. (|14|l . our approximate derived con- 
straint Eq. (|23|l for the red-shift range over which a structure 



of N particles may be resolved, can be rewritten as 



10" 



■N 



1/3 



(27) 



Thus, for example, at 



10 resolution constraints be- 



come relevant for all objects of less than N ~ 10 formed at 
or before a redshift of ten. 

Comparing this with Eq. (|25|l . in which the factor K,{e) 
has a priori very weak (logarithmic) dependence on N, leads 
to the conclusion that, for given e/i, two body coUisionality 
will be the dominant discreteness effect at small N, while for 
larger N it is the resolution limit associated with e which is 
the relevant one. Thus, for an e/i which remains fixed in co- 
moving coordinates, there is an inevitable trade-off between 
resolution and the introduction of spurious discreteness ef- 
fects. We have considered solely simulations with fixed co- 
moving softening, but the test can be applied without mod- 
ification (as described in detail in that last section) to any 
other kind of code. Indeed it could provide very useful con- 
straints on optimal smoothing strategies in codes with an 
adaptative smoothing, which indeed aim to optimize spa- 
tial resolution while keeping two b ody coUisionality under 
control (see e.g. iKnebe et all |20oO)). 

We finally remark that while the test discussed here 
can, as we have shown, provide constraints on the parame- 
ters which must be respected in order to reproduce the de- 
sired continuum limit, it does not show that satisfaction of 
these constraints guarantee the same property, i.e., the test 
provides necessary, but not sufficient, conditions to guar- 
antee the correctness of the numerical results. Firstly the 
test only constrains the choice of parameters in the strongly 
non-linear regime. This means, notably, that it cannot say 
anything about the appropriateness of the use of an e < £, 
which has been shown to introduce discreteness effects in 
the early time evolution, and which may or may not distort 
the sub sequent evolut ion (Splinte r et al. ([l998l): Jovce et al.l 
((20di); iRomeo et al] ((2008); Kncbe et all (|200oi ')'). Nor can 
it, as noted, provide const raints on the choic e of the start- 
ing redshift of simulations (|jovce et al.l (|2008l ) ; iKnebe et al.l 
(,200^)). Further, even in the non-linear regime, it is clear 
that considerable caution should be adapted in supposing 
the constraints derived from this test guarantee a faithful 
representation of the coUisionless evolution: it shows that 
halos with less than some number of particles will necessar- 
ily suffer from effects of discreteness which modify strongly 
their density profiles: as we have seen in our study, the use 
of a smoothing which is slightly too small leads at the end of 
the simulation to a completely incorrect profile. Further the 
characteristic physical scale of this modification is the size of 
the structure, which has no direct relation to the smoothing 
scale e. In the case of an inappropriately large e (simulation 
s7 above) we saw that a halo can even be very considerably 
larger than it should be. How the presence of such spurious 
clustering modifies the evolution of the whole system is very 
unclear. Indeed given that clustering in currently favored 
cosmological models is hierarchical in nature, the possibility 
that such error may feed through different scales is a major 
concern. 

We thank Frangois Sicard for useful discussions, and the 
two anonymous referees for many invaluable criticisms and 
suggestions. 
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APPENDIX A: FURTHER DETAILS ON 
SIMULATIONS 

Besides the choice of e on which we have focussed above, 
simulation with GADGET requires one to fix several other 
numerical parameters. In our choices we have followed the 
guidelines given by the GADGET user's guide, and also per- 
formed several tests indicate that our results appear to be 
reasonably stable with respect to these choices. Specifically 
we have considered: 

• ErrTolIntAccuracy, a dimensionless parameter which 
controls the accuracy of the time-step criterion. The value 
suggested by the GADGET user's guide is 0.025, and we 
have done numerical tests in various cases reducing it by a 
factor ten without any detectable difference in the results 
reported here. 

• MaxSizeTimestep, which specifies the maximum al- 
lowed time-step for cosmological simulations as a fraction of 
the current Hubble time. According to the GADGET user's 
guide a value of 0.025 is usually accurate enough for most 
cosmological runs. We have found our results to be stable in 
tests using considerably smaller values, down to as small as 
10~^ in some cases. 

• ErrTolTheta is the accuracy criterion (the opening an- 
gle 0) of the tree algorithm if the standard Barnes & Hut 
(BH) opening criterion is used. The suggested value is 0.7 
and we have considered values down to 0.1, again finding 
stable results in the tested cases. 

We have also performed some tests comparing different 
realizations of the initial conditions in various cases, again 
with stable results. 

We emphasize that despite these numerous tests, we 
have come to the conclusion through our analysis using the 
test studied in the paper that the four simulations with 
smaller e are in fact not numerically converged, and are 
characterized by very significant violations of energy con- 
servation due to poorly integrated two body collisions. Thus 
sufficient further extrapolation of the numerical parameters 
beyond what we have considered must lead to very different 
results in these cases. 

For completeness we give in Table lAll 1 the values of 
the numerical parameters used in the simulations reported 
in the body of the article. For the non-expanding simula- 
tion in open bound ary conditions we have used (see also 
ISvlos Labinil (|2012l ) for further details), a force softening 
£ = 0.007-R. This means that e is very much smaller at all 
times than the size of the structure. In addition we have 
chosen r; = 0.01, and have used the new GADGET cell 
opening criterion with a high force accuracy of ap = 0.001. 
E nergy is conserved to within less than 10^'^ up to 12rsc. 
In I Joyce et al] (|2009l ') extensive tests of the dependence on 
e were performed, for the much more constraining case of 
evolution from the same initial condition but with initial 
velocities set to zero. These tests lead to the conclusion that 
no sensible dependence is observed up to this time (notably 
in the density profile) unless the ratio e/Ra becomes large. 
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Simulation 


ETIA 


MaxTimeStep 


ErrTolTheta 


si 


0.025 


0.01 


0.7 


s2 


0.0025 


0.01 


0.7 


s3 
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0.01 


0.7 


s4 
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0.7 


s5 
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0.7 


s6 
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0.00001 


0.1 


s7 
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0.0001 


0.7 


slN5e4 


0.025 


0.01 


0.7 


slNleS 


0.025 


0.01 


0.7 


s5Nle3 


0.025 


0.0001 


0.7 


s6Nle3 


0.025 


0.00001 


0.1 



Table Al. Values of the three numerical parameters (ErrTolIn- 
tAecuraey, MaxTimeStep and ErrTolTheta) adopted in each of 
the simulations reported in the text. 
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